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ABSTRACT 

We study the dynamical stability against bar-mode deformation of rapidly and differentially rotating 
stars in the first post-Newtonian approximation of general relativity. We vary the compaction of the star 
M/R (where M is the gravitational mass and R the equatorial circumferential radius) between 0.01 and 
0.05 to isolate the influence of relativistic gravitation on the instability. For compactions in this moderate 
range, the critical value of j3 = T/W for the onset of the dynamical instability (where T is the rotational 
kinetic energy and W the gravitational binding energy) slightly decreases from ~ 0.26 to ~ 0.25 with 
increasing compaction for our choice of the differential rotational law. Combined with our earlier find- 
ings based on simulations in full general relativity for stars with higher compaction, we conclude that 
relativistic gravitation enhances the dynamical bar-mode instability, i.e. the onset of instability sets in 
for smaller values of f3 in relativistic gravity than in Newtonian gravity. We also find that once a triaxial 
structure forms after the bar-mode perturbation saturates in dynamically unstable stars, the triaxial 
shape is maintained, at least for several rotational periods. To check the reliability of our numerical 
integrations, we verify that the general relativistic Kelvin-Helmholtz circulation is well-conserved, in ad- 
dition to rest- mass energy, total mass-energy, linear and angular momentum. Conservation of circulation 
indicates that our code is not seriously affected by numerical viscosity. We determine the amplitude and 
frequency of the quasi-periodic gravitational waves emitted during the bar formation process using the 
quadrupolc formula. 

Subject headings: Gravitation — hydrodynamics — instabilities — relativity — stars: neutron — stars: 
rotation 



1. INTRODUCTION 

Stars in nature are usually rotating and subject to non- 
axisymmetric rotational instabilities. An exact treatment 
of these instabilities exists only for incompressible equi- 
librium fluids in Newtonian gravity, e.g. (Chandrasekhar 
1969; Tassoul 1978; Shapiro and Teukolsky 1983). For 
these configurations, global rotational instabilities arise 
from non-radial toroidal modes e %mip (m = ±1,±2,...) 
when j3 = T/W exceeds a certain critical value. Here if 
is the azimuthal coordinate and T and W are the rota- 
tional kinetic and gravitational potential binding energies. 
In the following we will focus on the m — ±2 bar-mode, 
since it is the fastest growing mode when the rotation is 
sufficiently rapid. 

There exist two different mechanisms and correspond- 
ing timescales for bar-mode instabilities. Uniformly rotat- 
ing, incompressible stars in Newtonian theory are secularly 
unstable to bar- mode formation when f3 > /3 SCC — 0.14. 
However, this instability can only grow in the presence of 
some dissipative mechanism, like viscosity or gravitational 
radiation, and the growth time is determined by the dis- 
sipative timescale, which is usually much longer than the 
dynamical timescale of the system. By contrast, a dy- 
namical instability to bar-mode formation sets in when 
P > /3dyn — 0.27. This instability is independent of any 
dissipative mechanisms, and the growth time is the hydro- 
dynamic timescale of the system. 

The secular instability in compressible stars, both uni- 



formly and differentially rotating, has been analyzed nu- 
merically within linear perturbation theory, by means 
of a variational principle and trial functions, by solving 
the eigenvalue problem, or by other approximate means. 
These techniques have been applied not only in New- 
tonian theory (Lynden-Bell and Ostriker 1967; Ostriker 
and Bodenheimer 1973; Friedman and Schutz 1978; Ipser 
and Lindblom 1990) but also in a post-Newtonian (PN) 
framework (Cutler and Lindblom 1992; Shapiro and Zane 
1998 for incompressible stars) and in full general relativity 
(Bonazzola, Frieben, and Gourgoulhon 1996, 1998; Ster- 
gioulas and Friedman 1998; Yoshida and Eriguchi 1999). 
For relativistic stars, the critical value of p sec depends on 
the compaction M/R of the star (where M is the gravita- 
tional mass and R the circumferential radius at the equa- 
tor), on the rotational law and on the dissipative mech- 
anism. The gravitational-radiation driven instability oc- 
curs for smaller rotation rates, i.e., for values (5 SCC < 0.14, 
in general relativity. For extremely compact stars (Ster- 
gioulas and Friedman 1998; Yoshida and Eriguchi 1999) 
or strongly differentially rotating stars (Imamura et al. 
1995), the critical value can be as small as /3 SCC < 0.1. By 
contrast, viscosity drives the instability to higher rotation 
rates (3 SCC > 0.14 as the configurations become more com- 
pact (Bonazzola, Frieben, and Gourgoulhon 1996, 1998; 
Shapiro and Zane 1998). 

Determining the onset of the dynamical bar-mode in- 
stability, as well as the subsequent evolution of an unsta- 



1 Department of Earth and Space Science, Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan 

2 Fortner Fellow 

3 Department of Astronomy and NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801 



2 



SAIJO, SHIBATA, BAUMGARTE & SHAPIRO 



ble star, requires a fully nonlinear hydro-dynamic simu- 
lation. Simulations performed in Newtonian theory, e.g. 
(Tohline, Durisen, and McCollough 1985; Durisen et al. 
1986; Williams and Tohline 1988; Houser, Centrella, and 
Smith 1994; Smith, Houser, and Centrella 1995; Houser 
and Centrella 1996; Pickett, Durisen, and Davis 1996; 
Toman et al. 1998; New, Centrella, and Tohline 2000) have 
shown that /3dyn depends only very weakly on the stiffness 
of the equation of state. Once a bar has developed, the 
formation of spiral arms plays an important role in redis- 
tributing the angular momentum and forming a core-halo 
structure. Recently, it has been shown that, similar to the 
onset of secular instability, /3dyn can be smaller for stars 
with a higher degree of differential rotation (Tohline and 
Hachisu 1990; Pickett, Durisen, and Davis 1996). 

To date, the dynamical bar-mode instability has been 
analyzed mainly in Newtonian theory, since until quite 
recently a stable numerical code capable of performing re- 
liable hydrodynamic simulations in three dimensions plus 
time in full general relativity has not existed. Some recent 
developments, however, have advanced the field signifi- 
cantly. New formulations of the Einstein equation based 
on modifications of the standard (3+1) system of equations 
have resulted in codes which have proven to be remarkably 
stable over many dynamical timescales, e.g. (Shibata and 
Nakamura 1995; Baumgarte and Shapiro 1999). In addi- 
tion, gauge conditions which allow long-time stable evolu- 
tion for rotating, self-gravitating systems and are manage- 
able computationally have been developed, e.g. (Shibata 
1999b). 

The purpose of this paper is twofold. We verify that the 
point of onset of dynamical bar mode formation, as mea- 
sured by /?, decreases with increasing compaction, and we 
furthermore show that in unstable configurations, the bar 
persists for at least several rotational periods. 

In a previous paper Shibata, Baumgarte, and Shapiro 
(2000) performed simulations of rapidly and differentially 
rotating neutron stars in full general relativity. They em- 
ployed the relativistic code of Shibata (1999a) to study the 
onset and growth of the dynamical bar-mode instability in 
relativistic stars. They found that for compact stars with 
M/R > 0.1, the onset of dynamical instability occurs at 
fldyn ~ 0.24 — 0.25, somewhat smaller than the Newto- 
nian value of /?d yn ~ 0.26. In principle, this reduction of 
/3dyn could be caused by effects of either general relativity 
or differential rotation. In order to isolate the two effects 
one could study sequences of varying compaction Mj R or 
sequences of varying degree of differential rotation. We 
choose the former in this paper and focus on the transi- 
tion to relativistic gravitation as the compaction increases 
to moderate values. 

In fully relativistic evolution codes the Courant condi- 
tion for the gravitational fields limits the size of the nu- 
merical timesteps due to the speed of light 4 . As a conse- 
quence, the ratio between the dynamical timescale for the 
matter and the Courant timestep for the fields scales ap- 
proximately as {M / R)^ 1 / 2 , which makes the calculation 
prohibitively slow for small compactions. In this paper, 
we therefore adopt the first order post-Newtonian (1PN) 
code of Shibata, Baumgarte, and Shapiro (1998) to study 



the effect of general relativity on /3dyn for small and mod- 
erate compactions. In this formalism, the relativistic evo- 
lution equations for the gravitational fields reduce to el- 
liptic equations, which are solved on individual timesliccs 
and hence have no Courant condition. The numerical 
timesteps are now limited by the Courant condition for 
the hydrodynamical evolution equations alone and there- 
fore scale with the dynamical timescale. We evolve models 
with compactions M/R between 0.01 and 0.05 for several 
rotational periods to determine their stability. In cases in 
which a bar forms we also follow the nonlinear growth and 
saturation of the instability. 

Should the star retain its bar-like shape for many rota- 
tional periods, then it would emit a quasi-periodic gravita- 
tional wave signal. Such a radiator would be potentially a 
very promising source for the gravitational wave interfer- 
ometers currently under construction. This prospect has 
prompted many researchers to study the evolution of the 
bar-mode instability, most recently New, Centrella, and 
Tohline (2000) and Brown (2000).' They have performed 
long-duration Newtonian simulations of dynamically un- 
stable stars that form a bar and then spiral arms. These 
configurations eject small amounts of mass and then settle 
down to triaxial stars. For the early stage of the evolu- 
tion, the results of these two groups are very similar, but 
for later times some differences arise. According to New, 
Centrella, and Tohline (2000), the bar gradually decays 
after a few rotational periods, probably due to numerical 
errors associated with the unphysical motion of the cen- 
ter of mass, whereas Brown (2000) reports that the star 
remains bar-like for more than 10 rotational periods, and 
presumably until gravitational radiation ultimately drives 
the decay of the bar. 

Given these discrepancies, which are probably due to 
numerical inaccuracies, it seems desirable to develop ad- 
ditional criteria to verify the late-time reliability of nu- 
merical codes. In particular, it is possible that numerical 
viscosity associated with finite-differencing could lead to 
an artificial decay of bar-modes. In this paper we evalu- 
ate the conservation of circulation, which is violated in the 
presence of viscosity (numerical or otherwise) . We present 
a method for computing the circulation in general rela- 
tivity and demonstrate that, in our numerical code, the 
circulation is well conserved. 

This paper is organized as follows. In Sec. 2 we present 
the basic equations of our 1PN formulation of general rel- 
ativity and describe our initial data in Sec. 3. We discuss 
our numerical results in Sec. 4, focusing on the dynami- 
cal stability of differentially rotating stars. In Sec. 5 we 
demonstrate the conservation of circulation in our sim- 
ulations and discuss the emitted gravitational wave sig- 
nal in Sec. 6. In Sec. 7 we briefly summarize our find- 
ings. Throughout this paper, we use the geometrized units 
with G = c = 1 and adopt Cartesian coordinates (x,y, z) 
with the coordinate time t. Greek and Latin indices take 
(t, x, y, z) and (x, y, z), respectively. 

2. BASIC EQUATIONS 

In this section, we briefly review the 1PN formalism of 
Shibata, Baumgarte, and Shapiro (1998). We solve the 

4 This limitation could in principle be avoided by using an implicit finite difference scheme. However, such schemes are more complicated than 
explicit schemes, and have not yet been implemented for fully relativistic hydrodynamics in three spatial dimensions 
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fully relativistic equations of hydrodynamics, but neglect 
some higher-order PN terms in the Einstein field equa- 
tions. 

2.1. The field equations 

Define the spatial projection tensor = g^ v + n^n v , 
where g^ is the spacetime metric, n M = (1/a, —(3 % /a) the 
unit normal to a spatial hypersurface, and where a and 
/?* are the lapse function and shift vector. Within a 1PN 
approximation, the spatial metric gij = 7^ may always be 
chosen to be conformally flat 

7y=-0 4 %, (1) 
where ip is the conformal factor (see Chandrasekhar 1965; 
Blanchet, Damour, and Schafer 1989). The spacetime line 
element then reduces to 

ds 2 = (-a 2 + P k /3 k )dt 2 + 2(3idx l dt 

+tP 4 5 ij dx i dx j . (2) 

We adopt maximal slicing, for which the trace of the ex- 
trinsic curvature Kij vanishes, 

K = ^Kij = 0. (3) 

The 1PN field equations for the five unknowns tp, a and 
(3 l can then be derived conveniently from the (3+1) for- 
malism. 

Since the spatial metric is conformally flat, the trans- 
verse part of its time derivative vanishes. The transverse 
part of the evolution equation of the spatial metric there- 
fore relates the extrinsic curvature to the shift vector, 



2aip 



- 4 K- 



(4) 



Inserting Eq. (4) into the momentum constraint equation 
then yields an equation for the shift vector f3 l 



S a A(3 l + ^dAp 1 = WiraJ, + (dj In ^ 
x (dip* + SuS^dkP 1 - ^S{diP 



(5) 



where A = 5 lJ didj is the flat space Laplacian and Ji = 
—n^'j^T^ is the momentum density. In the definition of 
Ji, T^ v is the stress energy tensor. 

The conformal factor tp is determined from the Hamil- 
tonian constraint 

AtP = -2^p H -\^K l3 K l \ (6) 
o 

where pn = n^n^T^ is the mass-energy density measured 
by a normal observer. 

Maximal slicing implies dtK = 0, so that the trace of 
the evolution equation for the intrinsic curvature yields an 
equation for the lapse function a, 

7 



A(aV>) = 27raV 5 (Pff + 2S) + -aip K^K 10 , (7) 

8 

7jfcT jfc . We also use Eq. (6) to derive this 



where S 
equation. 

Equations (5) - (7) determine the fully nonlinear rel- 
ativistic metric for maximal slicing within the conformal 
flatness approximation. None of these equations involve 
time derivatives, so that in a numerical finite difference 
implementation the Courant condition is no longer cou- 
pled to the speed of light. Since the conformal flatness ap- 
proximation introduces errors at higher order than 1PN, 



it is reasonable to neglect terms in these equations which 
are also higher than 1PN order. In particular, we discard 
[djiHa/iP 6 )}}^ + 8 a V k d k l3 l - 26fd l l3 l /3) in Eq. (5), 
tp 5 K i:j K^/8 in Eq. (6), and 7atp b K tJ K 1 ^ /8 in Eq. (7). 
With these simplifications, the source terms on the left 
hand sides of Eqs. (5) - (7) become compact. As a con- 
sequence, we can solve these equations numerically very 
accurately with outer boundary conditions set at fairly 
small distances outside the matter, and hence on fairly 
small numerical grids. 

The equation for the shift, 



6 u A(3 l + ^didtP 1 = 16W 4 , 



(8) 



can be further simplified by introducing a vector Bi and a 
scalar \ according to 



AB t = inaJi, 
Ax = -AiraJiX 1 . 

The shift can then be computed from 

5 ij /3 j =4B i -hd iX + d i (B k x k )}, 



(9) 
(10) 

(11) 



and will automatically satisfy Eq. (8). The vector-type 
Poisson equation (Eq. (8)) for /3 Z has hence been reduced 
to four scalar-type Poisson equations for Bi and 

To summarize, we have reduced Einstein equations in a 
1PN formalism to six elliptic equations for the six variables 
(atp, tp, B t , x), 

A(m/;) = 27ro^W + 2S) = ^S al p, (12) 

Aip = -2irip 5 pH = 4^6"^,, (13) 

ABi = 47raJ 4 , (14) 

Ax = -AiraJ l x\ (15) 

These Poisson-type equations are solved imposing the fol- 
lowing boundary condition at outer boundaries 



aip = l--j S atf) d 3 x + 0(r 3 ), 

V = 1 — ^ / S^d 3 x + 0{r- 3 ), 
B x = --r / aJ x xd 3 x - ^- I aJ x yd 3 x 



+0(r- 4 ), 



B n 



J a,J y xd 3 x — ^ J a,J y yd 3 x 



B,= 



+0(r- 4 ), 

~ J aJ z zd 3 x + 0(r- 4 ), 
X= IJ a ^ l d 3 x + 0{r- 3 ). 



(16) 
(17) 

(18) 

(19) 
(20) 

(21) 



2.2. The matter equations 

For a perfect fluid, the energy momentum tensor takes 
the form 



r=/)(l|£|-|tiV|Pf, (22) 
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where p is the comoving rest-mass density, e the specific 
internal energy, P the pressure, and the four-velocity. 
We adopt a T-law equation of state in the form 

P = (Y - l) P e, (23) 

where Y is the adiabatic index which we set to be 2 in this 
paper. 

In the absence of thermal dissipation, Eq. (23), together 
with the first law of thermodynamics, implies a polytropic 
equation of state 

P = np 1+1/n , (24) 

where n = l/(Y— 1) is the polytropic index and k is a con- 
stant. Thus, if the matter satisfies a polytropic equation 
of state initially, the polytropic form of the equation of 
state is preserved during the subsequent evolution in the 
absence of shocks. 

From V^T 7 *" = together with the equation of state 
(Eq. (23)), we can derive the energy and Euler equations 
according to 



~dt 



+ 



d(p*Uj) 
dt 



d(p*u i v : ') 



+ 



+p lf u j f3 J , 



= —aip 6 P t i — p+avfa.- 



+ 



2p*UkUk 



lp,i 



(25) 



(26) 



where 



v = 



p* 



(27) 
(28) 

(29) 
(30) 
(31) 



(pe) 1/r c 
dx l u l 

pa^tp 6 , 
„ -(l+re)w*, 
Ui = (1 + Ys)ui. 

Note that we treat the matter fully relativistically; the 
1PN approximation only enters through simplifications in 
the coupling to the gravitational fields. Note also that we 
do not need to include an artificial viscosity, since we do 
not encounter any shocks in the simulations in this pa- 
per. As a consequence we also do not need to solve the 
continuity equation 

dp* d(p*v l ) 
dt dx l 

since in the absence of shocks it is equivalent to Eq. (25). 

The gravitational mass M, e.g. (Bowen and York 1980), 
rest mass Mo, proper mass M p , angular momentum J, ki- 
netic energy T, and gravitational binding energy W of a 
rotating star can be computed from 



0, 



(32) 



M = -- 



2tt 



/ 



(p + pe + P)(<V) 2 
U, = / pify/^gtfx = j P* 



drx, 



M p = 



J = 



J pu t (l + e)y/=gd 3 x = J p»(l + e)d 3 x, 
J Tly/=g~d 3 x = J (xJ y - yJ x )iP 6 d 3 x, 



(33) 
(34) 
(35) 
(36) 



T = ^ I HT! 



gd 6 x 



= \ J n(xJ y - yJ x )^ 6 d 3 x, (37) 

W = M p + T - M. (38) 
We also define the quadrupole moments as 



,x i x j d 3 x, (39) 



and the nondimcnsional, scale-invariant ratio (3 = T/W , 
which is very useful to characterize the dynamical stability 
against bar-mode deformation. 

Since we use a polytropic equation of state, it is con- 
venient to rescalc all quantities with respect to n. Since 
K n/2 dimensions of length, we introduce the following 
nondimcnsional variables 



t = n- ,l / 2 t, 



>/2 



T, X = K 



V2, 



y = k 



-n/2 



y, 



z = K ~ n / 2 z, n = n n / 2 n, (40) 

M = K' n ^ 2 M, R = K~ n / 2 R, C = KT n l 2 C 

where t is the proper time, Q the angular velocity of the 
star, and C circulation (see Sec. 5 for a definition). We 
also define the central rotation period as P c = 2tt/£L . 
Henceforth, we adopt nondimcnsional quantities, but omit 
the bars for convenience (equivalently, we set k = 1). 

3. INITIAL DATA 

To prepare initial data, we construct axisymmetric ro- 
tating stars in equilibrium and slightly perturb them. The 
equilibrium configurations arc obtained by solving the 
equation of hydrostatic equilibrium together with the field 
equations for metric (i.e., Eqs. (12) - (15)). For station- 
ary solutions, the 1PN Euler equation can be integrated 
to yield the 1PN Bernoulli equation 

1 



ln(l + r £ ) + - ln[a 2 - ijj 4 m 2 (n + (3 v f\ 



+ 



J vfu^dQ. = 



const, 



(41) 



where w = ^ x 2 + y 2 and /3 V = (xf3 y - y(3 x )/w 2 . 

Following previous studies (Komatsu, Eriguchi, and 
Hachisu 1989a, b; Cook, Shapiro, and Tcukolsky 1992, 
1994; Bonazzola et al. 1993; Salgado, Bonazzola, Gour- 
goulhon, and Haensel 1994; Shibata, Baumgarte, and 
Shapiro 2000), we adopt the differential rotation law 



F(Q) = = A 2 (n - O), 



(42) 



where A is a constant with dimension of length and ilo 
is the angular velocity on the rotational axis. With this 
choice, the hydrostatic equation (Eq. (41)) can be inte- 
grated analytically. In the Newtonian limit ( it* — > 1 and 
u v — > zu 2 fl) the corresponding rotational profile reduces 
to 

n - A2Qo (43) 

This equation implies that A determines the length scale 
over which fi changes, so that a smaller value of A yields 
a larger degree of differential rotation. As in Shibata, 
Baumgarte, and Shapiro (2000), we choose A — r ei where 
r e is the equatorial coordinate radius of the star, which 
corresponds to a moderate degree of differential rotation. 
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This choice allows us to compare directly with the results 
of Shibata, Baumgarte, and Shapiro (2000) and to focus 
on the effects of general relativity. Note that only uni- 
formly rotating stars with sufficiently stiff equations of 
state (n < 0.5) become dynamically unstable to bar for- 
mation before reaching the mass shedding limit, so that 
the effects of general relativity on /3d yn have to be studied 
for differentially rotating stars for typical configurations. 

We prepare initial conditions for various values of the 
compaction M/R and rotation rate. The latter is conve- 
niently parameterized by the deformation r p /r ei where r p 
is the polar coordinate radius. Since we adopt a 1PN ap- 
proximation all results are correct only up to the linear or- 
der in GM/c 2 R (or (v/c) 2 where v is a typical speed), im- 
plying that we can expect reliable results only for moderate 
compactions M/R <C 0(1). We therefore restrict our anal- 
ysis to compactions in the range between 0.01 and 0.05. 
In Shibata, Baumgarte, and Shapiro (2000), we found that 
fully rclativistic stars with M/R ~ 0.10 — 0.15 and A = r e 
become dynamically unstable against bar-mode formation 
when (3 > 0.25. Taking this result as a guide, we prepare 
rotating stars with (3 ~ 0.25, for which the corresponding 
deformation r p /r e takes values ~ 0.25 — 0.33. We label 
different initial data models with a Roman number and a 
Latin letter, e.g., Model II (c), where the Roman number 
labels the compaction, and the Latin letter the deforma- 
tion, hence /3. We tabulate the physical parameters of our 
initial value models in Table 1. 

In order trigger bar-mode formation in unstable models, 
we slightly perturb the density of the equilibrium models 
according to 

(2 2 \ 

1 + 0.1 x X ~ 2 V j. (44) 

This perturbation affects only the I xx and I yy compo- 
nents of the quadrupole moment, which change by approx- 



imately 1%. This perturbation can therefore be considered 
linear initially. 

Both for the construction of initial data and their sub- 
sequent evolution, we assume planar symmetry across the 
equator, and solve the equations on a uniform grid of typ- 
ical size 101 x 101 x 51. In the axisymmctric initial config- 
uration, the star's major and minor axes are covered by 40 
and 10 — 13 grid points. Dynamically unstable stars with 
high values of [3 form bars and eject mass. To avoid signif- 
icant mass outflow across the outer boundaries during the 
simulation, we use a larger grid of 141 x 141 x 71 gridpoints 
for these models. We terminate any simulation when 1% 
of the total rest mass has been ejected from the numerical 
grid or the evolution time reaches around 8P C . Our longest 
runs took about 18000 time steps, and consumed about 80 
CPU hours on a VX/4R vector-parallel computer at the 
National Astronomical Observatory of Japan. 

4. DYNAMICAL STABILITY OF DIFFERENTIALLY 
ROTATING STARS 

We evaluate the stability of the perturbed rotating stars 
by monitoring the distortion parameter 



± XX "T J-yy 

which is a measure of the magnitude of the bar-mode per- 
turbation. In Fig. 1, we show r\ as a function of time for 
all our models. When the star is dynamically unstable, 
rj grows exponentially up to a saturation level at which 
rj = O(l). Once the perturbation has saturated, the max- 
imum value of rj remains nearly constant on dynamical 
timescales implying that the bar structure persists. For 
stable stars, the maximum value of rj <C 1 remains ap- 
proximately constant throughout the evolution. 

The early exponential growth in unstable stars can be 
seen more clearly in Fig. 2, where we plot \rj\ as a func- 



Table 1 

Differentially rotating stars in equilibrium. 



Model 


r p /r e 


Pmax 


Pc 


Pe 


M 


M 


J 


T/W 


R/M 


stability 


1(a) 


0.250 


0.0117 


37.40 


83.52 


0.113 


0.123 


0.0264 


0.265 


20.01 


unstable 


1(b) 


0.275 


0.0125 


36.03 


80.50 


0.110 


0.120 


0.0245 


0.259 


19.94 


unstable 


1(c) 


0.300 


0.0135 


34.97 


78.27 


0.107 


0.116 


0.0221 


0.249 


20.00 


unstable 


1(d) 


0.325 


0.0148 


33.96 


75.97 


0.109 


0.103 


0.0198 


0.238 


20.01 


stable 


II (a) 


0.250 


0.00655 


52.39 


111.6 


0.0711 


0.0746 


0.0128 


0.271 


34.34 


unstable 


II (b) 


0.275 


0.00753 


50.21 


107.1 


0.0700 


0.0735 


0.0120 


0.263 


34.04 


unstable 


II (c) 


0.300 


0.00821 


47.19 


100.9 


0.0709 


0.0746 


0.0116 


0.252 


32.37 


stable 


II (d) 


0.325 


0.00901 


46.19 


98.48 


0.0675 


0.0710 


0.0102 


0.240 


32.75 


stable 


III (a) 


0.250 


0.00359 


72.64 


150.6 


0.0417 


0.0428 


0.00566 


0.273 


61.36 


unstable 


III (b) 


0.275 


0.00381 


70.78 


146.7 


0.0401 


0.0412 


0.00511 


0.265 


62.46 


unstable 


III (c) 


0.300 


0.00430 


67.95 


140.4 


0.0393 


0.0404 


0.00470 


0.254 


61.73 


stable 


III (d) 


0.325 


0.00480 


66.15 


137.1 


0.0377 


0.0388 


0.00417 


0.241 


61.76 


stable 


IV (a) 


0.250 


0.00218 


94.58 


193.1 


0.0262 


0.0266 


0.00279 


0.275 


100.1 


unstable 


IV (b) 


0.275 


0.00236 


91.39 


186.6 


0.0256 


0.0260 


0.00258 


0.267 


100.2 


unstable 


IV (c) 


0.300 


0.00265 


88.21 


180.7 


0.0248 


0.0251 


0.00232 


0.255 


100.4 


stable 


IV (d) 


0.325 


0.00296 


85.81 


174.6 


0.0238 


0.0242 


0.00207 


0.242 


100.1 


stable 



Note. — Pmax: maximum rest-mass density; P c : rotational period along the rotational axis; P c : 
rotational period at the equator; M: gravitational mass; M : rest mass; J: angular momentum; T: 
rotational kinetic energy; W: gravitational binding energy; R: equatorial circumferential radius. 
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tion of time on a logarithmic scale. We can determine the 
growth time r g and the oscillation period r D by fitting to 
a function 

ry = r ?0 10 t/T9 cos(27ri/r + ^o), (46) 
where r/o and ip are constants. We tabulate T g and t d 
for the unstable stars in Table 2. Interestingly, t de- 
pends only very weakly on M/R and 0, which agrees with 
the fully relativistic findings of Shibata, Baumgarte, and 
Shapiro (2000). 

Table 2 

t and t g in the early stage of bar formation. 



Model 


To/Pc 


r a /Pc 


1(a) 


1.30 


1.54 


1(b) 


1.31 


2.17 


1(c) 


1.29 


3.52 


II (a) 


1.30 


1.57 


II (b) 


1.31 


2.30 


III (a) 


1.31 


1.60 


III (b) 


1.31 


2.45 


IV (a) 


1.30 


1.61 


IV (b) 


1.30 


2.63 



Note. — t 
bation. t : 



oscillation period of bar-mode pertur- 
growth time of bar-mode perturbation 



(see Eq. (46)). 

Stable stars may show signs of an exponential growth 
very early on, but their distortion parameter always re- 
mains very small i)<0(l). As a criterion for stability, we 
therefore check whether r\ follows an exponential growth as 
in Eq. (46) up to large values r\ = 0(1). Judging from this 
criterion, Models I (a), (b), (c), II (a), (b), III (a), (b), and 
IV (a) and (b) are dynamically unstable. We summarize 
these results in Fig. 3, where we denote stable and unsta- 
ble models in a M/R versus plane. The critical value 
/?dyn approaches ~ 0.26 in the Newtonian limit M/R — ► 0, 
and decreases to ~ 0.25 for M/R — 0.05. Our Newtonian 
value differs slightly from the result for uniformly rotating, 
incompressible stars (0 = 0.27), which is most likely due 
to the differential rotation in our models. This result con- 
firms the fully relativistic results of Shibata, Baumgarte, 
and Shapiro (2000), which we have also included in Fig. 3. 
There, we showed that 0d yn decreases to < 0.24 for large 
compactions, M/R > 0.1. 

In Figs. 4 - 7, we show contours of the density p* in the 
equatorial plane for the final stages of the dynamical evo- 
lution. These plots clearly exhibit a triaxial structure for 
the unstable models, while for stable models the density 
distribution hardly changes during the evolution. 

Once the bars in dynamically unstable models have sat- 
urated, they persist with very little change for several ro- 
tation periods (see, e.g., Models I (c) and IV (b) in Fig. 1). 
In realistic systems, gravitational radiation reaction would 
provide a dissipation mechanism which would cause a de- 
cay of the bars. Using the quadrupolc formula 



E ~ —£lo(I xx - Iyy) 



M 2 R 4 n 6 r 1 



(47) 



(see, e.g., Lai and Shapiro 1995), we can estimate the ra- 
tio between the radiation reaction timescale tqw an d the 



central rotation period as 

R 

M 



T GW 



1 



5/2 (M/R^ 3/2 



(48) 



Eq. (48) implies that as long as the star is not extremely 
compact M/R = 0(1), the timescale of the radiation re- 
action is much longer than the dynamical timescale of the 
system, even if its rotation rate is near break-up, fig ~ 
M/R 3 , and the star is highly deformed, r\ = 0(1). In our 
simulation we choose M/R < 0.05 so that tqw / Pc ^> 1. 
Gravitational radiation reaction, which is not present at 
1PN order, is therefore irrelevant for the calculations in 
this paper. 

5. CONSERVATION OF CIRCULATION 

Considerable effort recently has gone into determining 
whether or not a bar, once it has saturated, persists for 
many rotational periods, or whether it decays very soon 
(see, e.g., New, Centrella, and Tohline 2000; Brown 2000). 
Differences in the results are probably due to numerical 
errors, and possibly caused by the presence of numerical 
viscosity associated with the finite differencing. In this 
section we describe how the conservation of circulation in 
general relativity can be used to check for the presence 
of numerical viscosity and to establish the reliability of 
long-time simulations. 

According to the Kelvin-Hclmholtz theorem, the rela- 
tivistic circulation 



C(c) 



hu^X^da, 



(49) 



is conserved in isentropic flow along an arbitrary closed 
curve c (see Carter 1979; Landau and Lifshitz 1982). Here, 
h = 1 + e + P/pis the specific enthalpy, a is a Lagrange 
parameter which labels points on the curve c, and A M is 
the tangent vector to the curve c [i.e., A M = (d/da)^]. 
Conservation of C can be verified by computing 

^-C(c) = ldau v V v {hu„\ ti ) 

« T Jc 

= jda [\»u u \7 u {hu^) + {hu ll )u"V 



i 



= - j> daX^V^h 

= 0. (50) 

Here, to derive the third line from the second line, we use 
the fact that — (d/dr)^ and A^ are coordinate basis 
vectors, and thus commute according to 

(51) 

1 and the Euler equation 

« A V A (H) = -V„/i, (52) 
to obtain the fourth line. Note that it is the derivative of 
C(c) with respect to the proper time r that vanishes, so 
that the circulation has to be evaluated on hypersurfaces 
of constant proper time as opposed to constant coordinate 
time. 

We check the conservation of circulation for three cases, 
namely the unstable Models I (b) and (c), and the stable 



We also have used u^u^ 
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Model I (d). For Model I (b) we evaluate the circulation 
along three closed curves, and for Models I (c) and (d) 
for four. Each loop is located in the equatorial plane, and 
is initially aligned with a constant density contour. For 
Model I (b) we choose loops which intersect x/r e = 0.5, 
0.625, and 0.75 on the x axis, and for Models I (c) and (d) 
an additional loop which intersects x/r e =0.875 and y = 0. 
We follow the evolution of each loop with the help of La- 
grangian tracers, whose trajectories are computed from 



dx l 
~dt 



dr 

~dt 



(53) 



The number of test particles representing each loop is 
80 — 140 depending on its initial location. We use a first 
order numerical scheme to integrate Eqs. (53) forward in 
time, once the (Eulerian) hydrodynamic flow has been de- 
termined. We evaluate the circulation along a loop c at 
a proper time t by interpolating the hydrodynamic vari- 
ables in both space and time to the current location of the 
Lagrangian tracers. 

Fig. 8 shows that the circulation is well conserved in 
all three Models, indicating that numerical viscosity only 
has very small effects in our code. In Fig. 9, we also show 
the location of the loops (Lagrangian particles) at the final 
time steps. As expected, the curves for the dynamically 
unstable stars (Models I (b) and (c)) become deformed, 
while the curves for the dynamically stable star (Model 
I (d)) remain close to spherical. The outermost loop in 
Fig. 8 (a) seems to indicate a small violation of the con- 
servation of circulation. However, this loop is close to the 
surface of a star which forms spiral arms, and is hence 
strongly deformed (see Fig. 9). The representation of this 
loop with Lagrangian tracers was therefore not sufficient 
to accurately evaluate its circulation. We tabulate the rel- 
ative errors in the circulation in Table 3. Except for the 
outermost loop in Model 1(b), the circulation is conserved 
up to ~ 1% in our simulations. 

We would like to emphasize that circulation is conserved 
in relativity even in the presence of gravitational radiation 
(as long as the fluid flow is isentropic and no shocks form). 
As a consequence, conservation of circulation can be used 
as a very strong code test in fully relativistic simulations 



as well as Newtonian or post-Newtonian simulations which 
include gravitational radiation reaction terms. 

6. GRAVITATIONAL WAVES 

We compute approximate gravitational waveforms by 
evaluating the quadrupole formula, neglecting all PN cor- 
rections. In the radiation zone, gravitational waves can 
be described by a transverse-traceless, perturbed metric 
hf^ with respect to a flat spacetime. In the quadrupole 
formula, hj^ can be expressed as (Misncr, Thornc, and 
Wheeler 1973) 



■TT 



2 d 2 



rTT 



1 



TT 
ij 1 kk 



Sin I, 



(54) 



rdt 2 r y 3' 

where r is the distance to the source, and TT denotes the 
transverse-traceless projection 

1 

2 JyJ 
with 

P t 3 = 5/ - h.v?, n l = x l /r. (56) 



jTT p o p b j _J_p pab j 

1 ij r i r j 1 ab c\ ij 1 abi 



Choosing the direction of the wave propagation to be along 
the z axis, the two polarization modes of gravitational 
waves can be determined from 



(57) 



Note that this quantity contains second time derivatives 
of Iij T , which are difficult to evaluate numerically. We 
therefore rewrite Eq. (54) as 



h 



TT 



2 d i pT 



dt\ 12 



1 



" x J 



rTT 



ij ± kk 



(58) 



and use the continuity equation (Eq. (32)) to eliminate 
the time derivatives in lJl T 



It 7 



(p*v l x j + p*x l v j )d 3 x, 



(59) 



e.g. (Finn 1989). We are then left with having to compute 
only a first time derivative numerically, which can be done 



Table 3 

Relative error of the circulation. 



Model 


r/r e 


Initial value 


Final value 


Relative error 


1(b) 


0.500 


1.198 


1.203 


0.4% 


1(b) 


0.625 


1.651 


1.632 


1.1% 


1(b) 


0.750 


2.053 


1.899 


7.5% 


1(c) 


0.500 


1.173 


1.187 


1.2% 


1(c) 


0.625 


1.603 


1.600 


0.2% 


1(c) 


0.750 


1.989 


1.992 


0.2% 


1(c) 


0.875 


2.330 


2.356 


1.1% 


1(d) 


0.500 


1.129 


1.142 


1.2% 


1(d) 


0.625 


1.535 


1.536 


0.1% 


1(d) 


0.750 


1.903 


1.910 


0.4% 


1(d) 


0.875 


2.228 


2.257 


1.3% 



Note. — Conservation of circulation for Models I (a), (c) and (d). 
We set each of the final circulation value at the end point in Fig. 8. 
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much more accurately. For observers along the z-axis, we 
find 

r h+ Id.- ■ . . , 

-w-mrt 1 -- 1 ^ (60) 



In Fig. 10, we plot waveforms for Models I (b), (c), and 
(d). For the unstable Models I (b) and (c), we find, as 
expected, a quasi-periodic oscillation with growing ampli- 
tude during the early bar formation. For the stable Model 
I (d) , we find a periodic waveform with approximately con- 
stant amplitude. 

In all three models, the frequency of the periodic oscil- 
lations is approximately 

* *77QHz(i^Y (62) 
1.3P C \ Pc J 

Once a bar forms and reaches saturation in unstable stars, 
the oscillation period reduces to slightly smaller values 
(~ 1.18P C for Model I (b) and - 1.27P C for Model I (c)), 
and accordingly the frequency shifts to slightly larger val- 
ues. This shift in the frequency is caused by the significant 
deformation of the stars by the bars. 

The gravitational wave amplitude in Model I (b) is ap- 
proximately 

<— (^)(^)- 

(63) 

Both the frequency and gravitational wave amplitude are 
consistent with the fully rclativistic results in Shibata, 
Baumgarte, and Shapiro (2000). 

A detection of these signals by kilometer size laser- 
intcrferometric gravitational wave detectors like LIGO, 
e.g. (Thorne 1995) might be feasible, because the wave 
form is quasi-periodic, which significantly increases its ef- 
fective amplitude (see Lai and Shapiro 1995). Also, radia- 
tion reaction may gradually shift the frequency to smaller 
values, and hence into a regime where LIGO is more sen- 
sitive. The effect of radiation reaction on the evolution of 



bar modes is not very well understood, though, except in 
incompressible stars, e.g. (Chandrasekhar 1970; Lai and 
Shapiro 1995) and should be the subject of future studies. 

7. SUMMARY 

We perform 1PN simulations of rapidly and differen- 
tially rotating stars to investigate general rclativistic ef- 
fects on the dynamical bar-mode instability for small com- 
pactions M/R < 0.05. By combining these PN results 
with the fully relativistic simulations of Shibata, Baum- 
garte, and Shapiro (2000) for configurations of higher com- 
paction, we conclude that the critical value of (3 = /3dyn 
decreases with increasing M/R. Thus, relativistic gravita- 
tion enhances the bar-mode instability. 

We also describe how conservation of circulation can 
be used to check by how much a code is affected by nu- 
merical visocity. In the presence of significant numerical 
viscosity, long-time evolution calculations become very un- 
reliable and may lead, for example, to erroneous evolution 
of a saturated bar. We show that in our calculations the 
circulation is well conserved, implying that our code is 
at most very weakly affected by numerical viscosity. We 
present a method for computing the circulation which can 
be applied in Newtonian, post-Newtonian and fully rela- 
tivistic calculations. 
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5689), NSF Grants AST 96-18524 and PHY 99-02833 
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Abroad) and the hospitality of the Department of Physics 
at UIUC; T. W. Baumgarte is pleased to acknowledge sup- 
port through a Fortner Fellowship. Numerical computa- 
tions were performed on the VX/4R machines in the As- 
tronomical Data Analysis Center of National Astronomical 
Observatory of Japan, and on the FUJITSU- VX/S vector 
computer at Media Network Center, Waseda University. 
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Fig. 1. — Deformation parameter r\ as a function of t/P c for our sixteen different models (see Table I). Solid, dashed, dash-dotted and 
dotted lines denote Models (a), (b), (c) and (d), respectively. We terminate the simulations after t ~ 8P C , or when 1% of the total rest mass 
has escaped from the computational grid (Models I (a), (b), II (a) (b), III (a), and IV (a)). 

Fig. 2. — Same as Fig. 1, but we show \r]\ on a logarithmic scale for t/P c < 4. 

Fig. 3. — Summary of our dynamical stability analysis. All our models are plotted in a (3 versus M/R plane, with stable stars denoted by a 
circle and unstable stars by a triangle. The solid circles and triangles are the models studied in this paper; the open circles and triangles are 
the models explored in full GR by Shibata, Baumgartc, and Shapiro (2000). We conclude that the critical value of (3 = ftjyn slightly decreases 
with increasing compaction M/R. This trend is emphasized by the dotted line, which shows /3dyn as a function of M/R approximately. 

Fig. 4. — Final density contours for p„ in the equatorial plane for Models I. The contour lines denote densities p* = 1.3 i X 10~ 3 
(i = 1, . . . , 15) and at times (a) t/P c = 2.72, (b) t/P c = 3.66, (c) t/P c = 7.77, and (d) t/P c = 8.16. 

Fig. 5. — Same as Fig. 4 for Models II. The contour lines denote densities p» = 6.0 i X 10~ 4 (i = 1, ... ,15) at times (a) t/P c = 2.80, (b) 
t/Pc = 4.27, (c) t/Pc = 8.38, and (d) t/P c = 8.25. 

Fig. 6. — Same as Fig. 4 for Models III. The contour lines denote densities p* = 3.1 i X 10~ 4 (i = 1, . . . , 15), at times (a) t/P c = 2.91, (b) 
t/Pc = 4.25, (c) t/Pc = 8.04, and (d) t/P c = 7.93. 

Fig. 7. — Same as Fig. 4 for Models IV. The contour lines denote densities p* = 1.7 ix 10~ 4 (i = 1, . . . , 14), at times (a) t/P c = 2.77, (b) 
t/Pc = 7.72, (c) t/Pc = 9.86, and (d) t/P c = 9.73. 

Fig. 8. — Circulation as a function of proper time r for various loops in Models I (b), (c) and (d). The loops have an initial radius of 
r/r-c = 0.5 (solid line), 0.625 (dashed line), 0.75 (dash-dotted line) and 0.875 (dotted line). 

Fig. 9. — Location of Lagrangian test particles at the end of the simulations for Model I (b) (at t/P c = 7.37, 7.26 listed from the outer 
loop), Model I (c) (at t/P c = 11.06, 10.95, 10.84, 10.75) and Model I (d) (at t/P c = 7.67, 7.59, 7.51, 7.45). Each loop is plotted for a 
constant value of t, but different loops correspond to slightly different values of r. We also include density contours for p* = 10 — 4 p*max 
(solid lines), which nearly coincide with the surface of the star. The long dashed, dashed, dash-dotted, dotted lines denote the locations of 
the test particles that are initially located along loops which intersect x/r e = 0.875, 0.750, 0.625, and 0.500 on the x axis. 

Fig. 10. — Gravitational waveforms rh+/M (solid lines) and rh x /M (dashed lines) as seen by a distant observer located on the 2;-axis. 
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